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Abstract — We show how the standard (Stormer-Verlet) split- 
ting method for differential equations of Hamiltonian mechanics 
(with accuracy of order r 2 for a timestep of length r) can be 
improved in a systematic manner without using the composition 
method. We give the explicit expressions which increase the 
accuracy to order r 8 , and demonstrate that the method work 
on a simple anharmonic oscillator. 

Index Terms — Splitting-method, Hamilton-equations, Higher- 
order-accuracy, Symplecticity 



I. Introduction 

I HE Hamilton equations of motion constitute a system 
of ordinary first order differential equations, 
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where the ' denotes differentiation with respect to time t, 
and H = H(q,p). They can be viewed as the characteristic 
equations of the partial differential equation 
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p(q,p;t) = Cp(q,p;t), 



with L the first order differential operator 
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generating a flow on phase space. If H does not depend 
explicitly on t, a formal solution of (|2|l is 



p(q,p;t) =e tL p(q,p;0) 



(4) 



In most cases this expression remain just formal, but one 
may often split the Hamiltonian into two parts, H = Hi + 
H2, with a corresponding splitting C = C\ + £2 such that 
the flows generated by C\ and £ 2 separately are integrable. 
One may then use the Cambell-Baker-Hausdorff formula to 
approximate the flow generated by C. One obtains the Strang 
splitting formula (T], @ 



e 5T£ 2 e r£! e |T£ 2 _ e r£ + JjT 3 [2£i+£ 2 ,[£ 1 X 2 ]] + - 
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which shows that time stepping this expression with a 
timestep r provides an approximation with relative accuracy 
of order r 2 , exactly preserving the symplectic property of 
the flow. 
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This corresponds to the symplectic splitting scheme of 
iterating the process of solving 
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Here the last part of one iteration may be combined with the 
first part of the next, unless one deals with time dependent 
systems or wants to register the state of the system at the 
intermediate times. 

From a practical point of view the most interesting prop- 
erty of this formulation is that it can be interpreted directly in 
terms of physical processes. For instance, for Hamiltonians 
H(q,p) — T(p) + V(q), a standard splitting scheme is to 
choose Hi = T and H2 = V. In that case |6]l corresponds 
to a collection of freely streaming particles receiving kicks 
at regular time intervals r, these kicks being dependent of 
the positions q of the particles. I.e, we may think of the 
evolution as a collection of kicks and moves [3|. 

It is not clear that this is the best way to approximate or 
model the exact dynamics of the real system. For instance, 
why should the motion between kicks be the free streaming 
generated by T(p)l There are more ways to split the Hamil- 
tonian into two integrable parts Q; the best splitting is most 
likely the one which best mimics the physics of equation ([lj. 
Further, since this equation is not solved exactly by (|6| for 
any finite value of t we need not necessarily choose H2 to 
be exactly H — Hi as long as it approaches this quantity 
sufficiently fast as r —> 0. We will exploit this observation 
to improve the accuracy of the splitting scheme |6]l in a 
systematic manner. 

We are, of course, not the first trying to improve on 
the Stormer-Verlet splitting scheme. An accessible review 
of several earlier approaches can be found in reference (5). 
Neri (6 1 has provided the general idea to construct symplectic 
integrators for Hamiltonian systems. Forest and Ruth |7) 
discussed the explicit fouth order method for the integration 
of Hamiltonian equations for the simplest non-trivial case. 
Yoshida [8| worked out a symplectic integrator for any even 
order, and Suzuki [9| presented the idea of how recursive 
construction of successive approximants may be extended to 
other methods. 

II. Harmonic Oscillators 

For a simple illustration of our idea consider the Hamil- 
tonian 



H{p,q)= l -{p 2 



(7) 



whose exact evolution over a time interval r is 



cos t sin t 
- sin t cos t 
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Compare this with a kick-move-kick splitting scheme over the 
same time interval, with i/kick = \kq 2 and H mme — ^mp 2 , 
where k and m may depend on r. One full iteration gives 
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the exact evolution is reproduced. If we instead choose a 
move-kick-move splitting scheme, with H move = \fhp 2 and 
#kick = \kq 2 , one iteration gives 
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which becomes exact if we choose 

2 r T 
m = — tan — , k - 
t 2 



(11) 



(12) 



It should be clear that this idea works for systems of har- 
monic oscillators in general, i.e. for quadratic Hamiltonians 
of the form 



H( q ,p)= l -(p T Mp + q T Kq) 



(13) 



where M and K are symmetric matrices. For a choosen 
splitting scheme and step interval r there are always modified 
matrices M T = M + 0(t 2 ) and K T = K + C(r 2 ) which 
reproduces the exact time evolution. For systems where M 
and K are too large for exact diagonalization, but sparse, a 
systematic expansion of M T and K T in powers of r 2 could 
be an efficient way to improve the standard splitting schemes. 

III. Nonlinear systems 

For a more general treatment we consider Hamiltonians of 
the form 



H{q,p)= l -p T Mp + V{q) 



A series solution of the Hamilton equations in powers of r 
is 

<Ze = q° +p a T - \d a VT 2 - ld a (DV)T 3 + 0{T 4 ), 

2 o 

Vl =Pa- d a Vr - ^d a (DV)r 2 (15) 
+ d a {^DV- l -D 2 V^j t 3 + 0(t 4 ) 
Here we have introduced notation to shorten expressions, 
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D = P ad a 7 D = (d a V)d a , 



(16) 



where we employ the Einstein summation convention: An 
index which occur twice, once in lower position and once 
in upper position, are implicitly summed over all available 



values. I.e, M ab d D = i\^ b M ah db (we will generally use the 
matrix M to rise an index from lower to upper position). The 
corresponding result for the kick-move-kick splitting scheme 
is 
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As expected it differs from the exact result in the third order, 
but the difference can be corrected by introducing second 
order generators 
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-D 2 Vt 2 , 
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V 2 = —DVt 2 
24 



(18) 



to be used in respectively the move and kick steps. Special- 
ized to a one -dimensional system with potential V = \q 2 
this agrees with equation ( fT0] >. With this correction the kick- 
move-kick splitting scheme agrees with the exact solution to 
4 th order in r, but differ in the r 5 -terms. We may correct the 
difference by introducing fourth order generators, 
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(19) 
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Specialized to a one-dimensional system with potential V 



7jq this agrees with equation (lOi. With this correction 



the kick-move-kick splitting scheme agrees with the exact 
solution to 6 th order in r, but differ in the r 7 -terms. We may 
correct the difference by introducing sixth order generators, 
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where we have introduced 



A = (d a v)(d b v)(d c v)d a d b d c . 



(20) 
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(14) Specialized to a one-dimensional system with potential V = 



\q 2 this agrees with equation (II Ok. With this correction 



the kick-move -kick splitting scheme agrees with the exact 
solution to 8 th order in r, but differ in the r 9 -terms. One 
may continue the correction process, but this is probably 
well beyond the limit of practical use already. 

IV. Solving the move steps 

Addition of extra potential terms V —> V e s = V + V2 + 
V4 + . . . is in principle unproblematic for solution of the kick 
steps. The equations, 



q a = 0, Pa = -d a V sS {q), 



(22) 



can still be integrated exactly, preserving the symplectic 
structure. The situation is different for the kinectic term 
T — > T e ff = T+T-2+T4 + - ■ ■ , since it now leads to equations 

.„ d 
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T eS (q,p), Pa = -d a T e «{q,p), (23) 



which is no longer straightforward to integrate exactly. Al- 
though the problematic terms are small one should make sure 
that the move steps preserve the symplectic structure exactly. 
Let q,p denote the positions and momenta just before the 
move step, and Q, P the positions and momenta just after. 
We construct a generating function fl0|-[|12| G(q,P;r), 
with 
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This preserves the symplectic structure; we just have to 
construct G to represent the move step sufficiently accurately. 
Consider first the case without the correction terms. The 
choice G = q a P a + \P a P a r gives 



Q a = q a + P a T, Pa = P a , 



(25) 



and find the first terms in the expansion to be 

Go =q a P a , 

Gl = \P a Pa, 
G 2 = 0, 
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which is the correct relation. Now add the T2-term to the 
move step. To order r 4 the exact solution of equation (23 i 
becomes 



(26) 



Q a = q a +p a r- \d a DV r 3 - ^-d a D 2 V r\ 
24 

Pa = Pa + ^d a D 2 Vr 3 + ±d a D 3 Vr\ 



Compare this with the result of changing 



G^G- ^V 2 Vt 3 - ^V 3 Vt 4 , (27) 



where V = P a d a . The solution of equation (24 1 change from 
the relations d25l) to 



q a + P a r - -d a VVr 3 - -d a V 2 V r 4 , (28) 
6 8 



Pa = P a - ^d a V 2 VT 3 - ^d a V 3 Vr\ 



(29) 



Since V is linear in P, equation (29 1 constitute a system 
of third order algebraic equation which in general must be 
solved numerically. This should usually be a fast process for 
small r. An exact solution of this equation is required to 
preserve the symplectic structure, but this solution should 



also agree with the exact solution of ( 23 1 to order r 



This may be verified by perturbation expansion in r. A 



perturbative solution of equation ( 29 1 is 



— 63 DVDV 2 + 25 VDVDV) V. 

V. EXPLICT COMPUTATIONS 
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Fig. 1. This figure illustrate how well energy is conserved with the various 
splitting schemes. The qualities plotted is (H — |)/T ,n for r = 0.2 
(squares), r = 0.1 (triangles) and t = 0.05 (lines). Here m = 2 for the 
Stormer-Verlet scheme (dotted line), m = 4 for the T 2 -corrected generators 
(dash-dotted line), m = 6 for the r 4 -corrected generators (dashed line), 
and m = 8 for the T 6 -corrected generators (fulldrawn line). Each plotted 
quantity is essentially the value of the next correction at the visited point 
in phase space. Since the plot is taken over the last half of the 16 th period 
the figure also give some indication of how well the exact oscillation period 
is reproduced by the scheme. The deviation is quite large for the Stormer- 
Verlet scheme when r = 0.2; to avoid cluttering the figure we have not 
included these points. 
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It remains to demonstrate that our algorithms can be ap- 
plied to real examples. We have considered the Hamiltonian 
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(32) 



which inserted into ( 28 i reproduces the full solution ( 26 1 to 
order t 4 . 

This process can be systematically continued to higher is a nonlinear oscillation with H constant equal to |, and 
orders. We write the transformation function as 



G(r) = J2G n T n , 



(30) 



with initial condition q(0) = 0, p(Q) = 1. The exact motion 
is a no 
period 
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Here B(x,y) — T(x)T(y)/T(x + y) is the beta function. In 
figure ll] we plot the behaviour of (i? — |) /r 2+ ™ during the 
last half of the 16 th oscillation, for various values of r and 
corrected generators up to order t 6 (corresponding to n = 6). 



Energy preservation with the Stormer-Verlet splitting scheme 



Energy preservation with 7- 4 -corrected generators 





Fig. 4. This figure illustrate the long time behaviour for the r 4 -corrected 
scheme (through the last half of the 4 104 th period). Different timesteps r 
have an effect on the period of oscillation, but the preservation of energy 
remains stable for a very long time. 



Energy preservation with 7- 6 -corrected generators 



Fig. 2. This figure illustrate the long time behaviour (through the last half 
of the 16 th period) for the Stormer-Verlet scheme. Different timesteps r 
have an effect on the period of oscillation, but the preservation of energy 
remains stable for a very long time. 



Energy preservation with 7- 2 -corrected generators 





Fig. 5. This figure illustrate the long time behaviour for the r B -corrected 
scheme (through the last half of period 262 718). Different timesteps r have 
an effect on the period of oscillation, but the preservation of energy remains 
stable for a very long time (for high accuracy and very long runs the effect 
of numerical roundoff errors eventually becomes visible). 



Fig. 3. This figure illustrate the long time behaviour (through the first 
half of the 257* period) for the t 2 -corrected scheme. Different timesteps 
r have an effect on the period of oscillation, but the preservation of energy 
remains stable for a very long time. 

VI. Conclusion 

We have shown that it is possible to systematically 
improve the accuracy of the usual symplectic integration 
schemes for a rather general class of Hamilton equations. 
The process is quite simple for linear equations, where it 
may be useful for sparse systems. For general systems the 
method requires the solution of a set of nonlinear algebraic 
equations at each move step. To which extent an higher-order 
method is advantageous or not will depend on the system 
under analysis, and the wanted accuracy. As always with 
higher order methods the increased accuracy per step may 
be countered by the higher computational cost per step (TS) . 
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